Comparison of methods for estimating continuous distributions of relaxation times 
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The nonparametric estimation of the distribution of relaxation times approach is not as frequently 
used in the analysis of dispersed response of dielectric or conductive materials as are other immit- 
tance data analysis methods based on parametric curve fitting techniques. Nevertheless, such distri- 
butions can yield important information about the physical processes present in measured material. 
In this letter, we apply two quite different numerical inversion methods to estimate the distribution 
of relaxation times for glassy Lio.sLao.sTiOs dielectric frequency-response data at 225 K. Both 
methods yield unique distributions that agree very closely with the actual exact one accurately cal- 
culated from the corrected bulk-dispersion Kohlrausch model established independently by means 
of parametric data fit using the corrected modulus formalism method. The obtained distributions 
are also greatly superior to those estimated using approximate functions equations given in the 
literature. 
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Broadband dielectric (also known as immittance or 
impedance) spectroscopy is widely used to character- 
ize materials and to help understand the mechanisms 
involved in such challenging areas of condensed-matter 
physics as conductivity, molecular relaxation, liquid-glass 
transition etc. 1]. In this experimental technique an 
electrical property of the material is recorded as a func- 
tion of probing field frequency v. Data may be ex- 
pressed at one of the four specific immittance levels 
(j) the complex resistivity p(u>); (ii) the complex mod- 
ulus M(ui) = iuje$p{uj); {Hi) the complex permittiv- 
ity e = [Af(a;)] ; and (iv) the complex conductivity 
o{uj) = iLUEoe(u)) = [p(w)] _1 . Here, u is the angular 
frequency u> = 2ttv; Eq is the permittivity of free space; 
and i = y— T. 

Once a data set is acquired, it may be expressed 
at an appropriate immittance level and then analyzed 
to obtain valuable information about material pro- 
cesses. Often employed procedures that have been used 
to analyze frequency response data are (a) using the 
Kohlrausch- Williams- Watt (KWW) approached derived 
from stretched exponential behavior in the time do- 
main |2j, |3j ; (b) the Havriliak-Negami empirical expres- 
sion and (c) estimation of the distribution of relax- 
ation times (DRT) inherent in the data HHBHEH, 
an approach not as commonly employed as the other two. 
Unlike the KWW analysis of (a) , procedure (b) is a data 
fitting method that does not lead to added understanding 
of the physical processes presented in the experimental 
material. On the other hand, KWW analyses involve 
fitting models whose parameters are all of physical sig- 
nificance. Although they are useful for comparing fit pa- 
rameters for various materials at different state variable 
levels they are less appropriate for data involving several 
DRTs associated with different physical processes. 

The DRT approach of (c) is an elegant method for in- 
vestigating the contributions of relaxing units to the to- 



tal relaxation and for determining the influence of state 
variables on the relaxation. In the presence of differ- 
ent processes or broad relaxations, the DRT approach 
is superior to the parametric ones since (1) no a priori 
assumptions are needed, i.e., a sum of empirical expres- 
sions etc.; (2) the actual distributions in a given data 
set are initially unknown; (3) a DRT can be related to 
various physical parameters of the system; (4) and when 
there are two different overlapping relaxations present, 
their depencies on state variables would be easy to iden- 
tify and to observe the influence of the state variables 
on the distributions. As an example, the dynamic com- 
plexity of the relaxation system can be determined by 
estimating its DRT and thus establishing whether it is 
intrinsicly broadening (homogeneous) or a distribution 
of responses (heterogeneous) [To| . A distribution may be 
characterized as either discrete (composed of individual 
points) or continuous, and DRT analysis can unambigu- 
ously distinguish between these two possibilities 0,l3l- 
Recently, non-resonant spectral hole burning technique 
has been employed to resolve distinct continuous distri- 
butions experimentally in order to identify multiple re- 
laxing domains in materials [Tl|. 

In this letter, we compare the results of two different 
DRT inversion methods for analyzing a set of experimen- 
tal frequency-response data that involves a continuous 
distribution. We also compare the accuracy of two equa- 
tions for estimating appropriate distribution functions 
proposed by Bottcher and Bordewijk [5j. Although esti- 
mation of discrete-point distributions is not an ill-posed 
mathematical problem |?J , distribution estimation of con- 
tinuous distributions, the usual situation, is ill-posed. It 
is therefore particularly important to assess the utility 
and power of different DRT estimation procedures for a 
well-defined data situation. 



Experimental data, with M — 52 pc 
Lio.sLao.sTiOg (LLT) glass at 225 K [ij, 
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FIG. 1: (a) Real and (b) imaginary parts of the raw dielectric 
permittivity of LLT at 225 K. The solid line and solid points 
shows the full data, including electrode effects and that as- 
sociated with ohmic conductivity. The open-circle symbols 
show the predictions of a full complex nonlinear least-squares 
fit of this data set, as described in the text. The subscript 'R' 
in the axis labels indicates that the data are presented with- 
out any transformation (raw data). The inset shows the same 
data plotted at the complex resistivity level, and the arrow in 
the inset indicates the direction of increasing frequency. The 
vertical dashed line shows the division between the two types 
of response present and defines the critical radial frequency 

We. 



the complex resistivity and dielectric levels, were found 
to involve an appreciable component associated with elec- 
trode polarization effects. LLT conducts by ionic hopping 
and involves a finite dc resistivity, <tq = cr(Q). Further 
analysis of data for this material over a range of tem- 
peratures established that both ctq and the characteristic 
relaxation time of the dispersion of the bulk material, 
t , were thermally activated with To§ and t d having the 
same activation energies [l3| . 

Such behavior indicates that it is most appropriate to 
identify the bulk dispersive response of this material with 
a conductive-system dispersion of resistivity relaxation 
times, rather than with a dielectric-system distribution of 
permittivity relaxation times, one where <7o would be nat- 
urally interpreted as a leakage conductivity unrelated to 
the bulk dielectric dispersion process. Since conductive- 
system response has already been analyzed for this data 
set 0], and since it has been shown by data fitting 
that it may often be difficult to discriminate between fits 
of conductive-system and dielectric-system models when 
only a single data set is available[6|, we have elected to 
compare the two different DRT estimation procedures 
by determining their dielectric-system DRTs from the 
present data expressed at the complex permittivity level. 



FIG. 2: (a) Real and (b) imaginary parts of the dielectric 
permittivity of LLT at 225 K after elimination of the con- 
tributions of ohmic conductivity and electrode polarization 
effects and generation of virtually exact KD-model data us- 
ing LEVM. The axis identifiers are shown with a subscript 
'mce' denoting subtracted ohmic conductivity and electrode 
effects. 



The two analysis methods considered here will be des- 
ignated I and II. Method I involves a weighted nonlinear 
least squares approach for estimating dielectric distribu- 
tion strength points, g i; at corresponding relaxation-time 
values n, with 1 < i < N 7}. It allows either discrete or 
continuous DRTs to be estimated in terms of the {g,, t{\ 
values and their uncertainties, with the set of Tj's either 
taken fixed or free to vary. Better results are nearly al- 
ways obtained with Tj's taken free, as in the present work. 
In addition, the data may be in temporal response form 
or in the frequency domain involving complex response 
or either its real or imaginary part. An extensive fit- 
ting and inversion program named LEVM that includes 
Method I is available for free downloads^]. Method II 
is based on a constrained least-squares with the Monte 
Carlo procedure @- It leads to delta sequence distribu- 
tions 0] when applied to discrete DRTs @- Recently, 
a method rather similar to that of II has been indepen- 
dently proposed, one that uses nonparametric Bayesian 
statistics for solving similar inversion problems |16( . 

Since we arc interested in the dielectric DRT for the 
dispersive bulk relaxation process, it is important to elim- 
inate the contributions to the data arising from partly 
blocking electrode effects before estimating the DRT. To 
do so, a KWW response model, the KD, involving a 
strctchcd-cxponential shape parameter [3d , a characteris- 
tic relaxation time r Q , and a Ae strength parameter, was 
used for fitting the original full data with inclusion of 
free parameters to model the electrode effects, the high- 
frequency- limiting bulk dielectric permittivity, £oo, and 
do UJi- The fit was excellent and yielded the follow- 
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FIG. 3: Various estimates of unnormalized distribution of re- 
laxation times following from the data presented in Fig. [5] 
Here the implicit scaling factor for these distributions was 
Ae. The thick solid line is the unnormalized KD-model dis- 
tribution. The method I open-square symbols show points 
estimated using the real part of that data and the open-circle 
ones show those obtained from the imaginary part. The sizes 
of the symbols were determined by the estimated uncertain- 
ties of the fits. The many solid-stair lines show the results 
of the method II estimation procedure using the full com- 
plex data. The vertical dashed line here and in Fig. 4 shows 
the critical time constant. For comparison, the inset presents 
method II estimates using the raw data of Fig. Q 



ing rounded estimates for /3d, t , Ae, and e^: 0.547, 
2.63 U-s, 137, and 112, respectively. 

The precise fit values of these parameters were then 
used in LEVM, omitting those of a and the electrode 
ones, to generate a set of M = 300 data points repre- 
senting just the bulk part of the Kohlrausch response to 
eight significant figures or better. This data set is used 
below to estimate its DRT by the methods mentioned 
above. In addition, given only values of (3d, Ae, and a 
set of logarithmically distributed values of r for the range 
from about 0.96 ns to 10 ms, LEVM was employed to 
calculated highly accurate values for the KD DRT com- 
parison with the inversion estimate. 

The complex dielectric permittivity may be expressed 
in terms of a general DRT formalism, 



(e s - £oo) / 
J — ( 



g(lnr)dlnT 



1 



(1) 



where, Soo = e'(oo) and e s = e'(0) (the quantity Ae = 
£s — £oc is defined as the dielectric strength), and g(lnr) 
is the distribution function. For a delta sequence distri- 
bution Eq. Q leads to simple Debye response [T^ |. 
Both applied methods I and II are based on Eq. Q 
and are further described in Ref. Macdonald Q and 
Rcf. Tuncer ||, respectively. 

In Fig.^ the complex dielectric permittivity raw data 
for LLT are presented without transformation. As ev- 
ident in the inset of Fig^i the data include two dif- 
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FIG. 4: Comparison of several different normalized DRT es- 
timates. Identifications are the same as those in Fig. 3 for 
similar quantities. Shown are the exact KD-model DRT, the 
method I real (□) and imaginary (O) data estimates, the 
exponential smoothing of method II results (~k), and approx- 
imate DRT estimates calculated directly from the data of Fig. 
2 using Eq. (2) (- -) and Eq. (3) (- ■ -). The inset shows 
method II estimates using the raw data of Fig. compared 
with the scaled exact DRT (all divided by a factor of 100). 



ferent processes, with the right spur part representing 
low-frequency electrode polarization effects. The dashed 
vertical line (- -) in the inset indicates the approxi- 
mate crossover position (shown at 40 krads -1 ) from bulk 
dielectric system dispersion to conductivity and double- 
layer effects ^jj- Since all the open-circle fit points in 
the figure enclose their corresponding solid data points 
symmetrically, one may conclude that the fit is excellent. 

After we remove the contributions of the ohmic con- 
ductivity and electrode effects to the raw data, as de- 
scribed above, the pure dielectric-system dispersion is 
obtained and is presented in Fig.|2and denoted by e mce - 
This data set, implicitly involving the KD-model DRT, 
was next used to estimate the DRT by the inversion 
methods I and II. Some of these results are shown in 
Fig. and H The thick solid line is that of the KD 
DRT with fa = 0.54657. Note that the data of Fig. H 
contain neither systematic nor random errors and thus 
allow comparison of the utility of methods I and II with- 
out such confounding factors. It is striking that the two 
methods both yield very accurate estimates of the exact 
DRT. The precision of the estimates obtained by method 
I using the real part of the data (Fig.[3J) is the best of the 
results shown and is remarkably small, especially for the 
points at and to the left of the peak. We should also re- 
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member that increasing the number of randomly selected 
t values used in method II improves the DRT estimates; 
~ 25000 t values were used in the present work. 

Method II selects random r values over a wider range 
than those defined by the range of the original frequency 
window. The range of the original data (er) is about 
2 krads -1 < lu < 200 Mrads" 1 corresponding to 5 ns < 
t < 500 |xs, somewhat smaller than the r range following 
from the exact data of Fig. |2 as defined above. 

In order to illustrate the utility of method II, its DRT 
determined for the raw £r data is shown in the inset of 
Fig. 13 with solid vertical lines and is compared to the 
actual distribution. Note that the presence of electrode 
effects results in an added distribution with a peak at 
t = 100 [is. In addition, the distributions obtained from 
the £r and e mce data sets are nearly the same for r < 
10 |xs except for the presence of a small peak of the £r 
distribution estimate near r ~ 32 ns. This could possibly 
be due to the raw data where no a priori assumption 
is made of the presence of KD-model response (e mce ). 
Also note the effects of the relaxation-time cutoff for fast 
processes at t < 1 ps. 

To further emphasize the utility of the numerical in- 



version methods for estimating a DRT, we compare our 
results with those of Bottcher and Bordewijk y| in Fig. 21 
They derived approximate DRT expressions from the real 
and imaginary parts of the dielectric permittivity and 
their derivatives with respect to natural logarithm of an- 
gular frequency (lnaj). Two such approximation distri- 
bution functions g are listed below. 

gxQnc,- 1 ) = 2 £ "H( 7 rAe)- 1 , (2) 
gating 1 ) = -Ae^ 1 de' (lu)/ din lu, (3) 

where the Fig. Q data fit result for Ae, 136.93, is used 
along with e mce data set values. Clearly these expressions 
lead to broader distributions and to far less accurate DRT 
estimates than our inversion ones. In the inset of Fig. 01 
method II DRT estimates obtained from the raw data 
are again illustrated, together with those following from 
the application of Eqs. @ and 

In conclusion, two approaches for estimating DRT in 
conductive and dielectric systems are applied to experi- 
mental LLT dielectric permittivity data at 225 K. Both 
methods are capable of yielding well defined unique dis- 
tributions for a given data set. 
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